#ifndef PCPS_ACCUMULATOR_PROPERTIES_HEADER
#define PCPS_ACCUMULATOR_PROPERTIES_HEADER

#include <src/pcps.h>
#include <src/accumulator.h>
#include <src/operator.h>
#include <src/local_correlation.h>

namespace pcps {

  // forward declaration of thread accumulator
  template <class S, class PART, class REF> class PropertiesThreadAccum;

  //--------------------------------------------------------------------------------
  // pcps::PropertiesProcessAccum -- Accumulates correlation functions.
  //--------------------------------------------------------------------------------

  template <class S, class PART, class REF> class PropertiesProcessAccum : public pcps::ProcessAccum<S, PART, REF> {

    private:

      // data members
      const int _nsites;                        // number of sites in the lattice
      const int _configs_per_site;              // the number of configurations in a single site's Hilbert space
      std::vector<S> _process_sz;          // array to hold this process's contribution to < Sz_i > for each site i
      std::vector<S> _process_sz_dot_sz;   // array to hold this process's contribution to < Sz_i Sz_j > for each pair of sites i, j
      std::vector<S> _process_s_dot_s;     // array to hold this process's contribution to < S_i . S_j > for each pair of sites i, j
      std::vector<S> _process_ssss;        // array to hold this process's contribution to < S_i.S_j S_k.S_l>
      const pcps::action_value_object<S, pcps::OpRL<S>, PART> _avoRL;  // object for computing numerical coefficient
                                                                                  // for action of raising/lowering operator

    public:

      // give the corresponding thread accumulator access to private and protected members
      friend class pcps::PropertiesThreadAccum<S, PART, REF>;

      // constructor
      PropertiesProcessAccum(const pcps::Input & userinp, const pcps::Wavefunction<S> & wfn, const long int sample_length)
        : pcps::ProcessAccum<S, PART, REF>(userinp, sample_length),
          _nsites(userinp.nsites()),
          _configs_per_site(userinp.configs_per_site()),
          _avoRL(userinp.configs_per_site())
      {

        // initialize arrays
        const S z = pcps::zero<S>();
        _process_sz.assign(         size_t(_nsites)                         , z );
        _process_sz_dot_sz.assign(  size_t(_nsites*_nsites)                 , z );
        _process_s_dot_s.assign(    size_t(_nsites*_nsites)                 , z );
        _process_ssss.assign(       size_t(_nsites*_nsites*_nsites*_nsites) , z );

        this->reset();

      }

      // function to create a thread accumulator linked to this process accumulator
      pcps::ThreadAccum<S, PART, REF> * create_thread_accum(const int tid, pcps::Arena & arena);

      // function to reset the accumulator
      void reset() {
        pcps::xscal(                        _nsites, pcps::zero<S>(), &_process_sz.at(0),        1);
        pcps::xscal(                _nsites*_nsites, pcps::zero<S>(), &_process_sz_dot_sz.at(0), 1);
        pcps::xscal(                _nsites*_nsites, pcps::zero<S>(), &_process_s_dot_s.at(0),   1);
        pcps::xscal(_nsites*_nsites*_nsites*_nsites, pcps::zero<S>(), &_process_ssss.at(0),      1);
      }

      // function to print information during stochastic iterations
      void print_iter_info(const pcps::Input & userinp,
                           pcps::Wavefunction<S> & wfn,
                           const int iter,
                           std::vector<boost::shared_ptr<pcps::ProcessAccum<S, PART, REF> > > & accumulators) {

        if (this->myrank() == 0) {
          const std::string filename = this->get_file_name(iter);
          std::cout << boost::format("properties accumulator:     properties written to file \"%s\"") % filename << std::endl;
          this->write_to_file(filename);
        }

      }

      // function to get the filename to write the properties to
      std::string get_file_name(const int iter) const {
        std::string retval = (boost::format("%i") % iter).str();
        if (retval.size() < 8)
          retval.insert(retval.begin(), 8 - retval.size(), '0');
        retval = (boost::format("properties_iter_%s.txt") % retval).str();
        return retval;
      }

      // function to write properties to a file
      void write_to_file(const std::string & filename) const {

        // open the output file
        std::ofstream outfile(filename.c_str());
        if (!outfile.is_open())
          throw pcps::Exception( (boost::format("failed to open file \"%s\" for writing properties") % filename).str() );

        // write the Sz averages
        outfile << boost::format("%6s  %40s") % "site i" % "< Sz_i >" << std::endl;
        for (int i = 0; i < _nsites; i++) {
          outfile << boost::format("%6i  %40s")
                     % i
                     % pcps::format_number("%14.6f", _process_sz[i] )
                  << std::endl;
        }
        outfile << std::endl;

        // write the Sz_i Sz_j averages
        outfile << boost::format("%6s  %6s  %40s") % "site i" % "site j" % "< Sz_i Sz_j >" << std::endl;
        for (int i = 0; i < _nsites; i++)
        for (int j = 0; j < _nsites; j++) {
          outfile << boost::format("%6i  %6i  %40s")
                     % i
                     % j
                     % pcps::format_number("%14.6f", _process_sz_dot_sz[i*_nsites+j] )
                  << std::endl;
        }
        outfile << std::endl;

        // write the S_i.S_j averages
        outfile << boost::format("%6s  %6s  %40s") % "site i" % "site j" % "< S_i . S_j >" << std::endl;
        for (int i = 0; i < _nsites; i++)
        for (int j = 0; j < _nsites; j++) {
          outfile << boost::format("%6i  %6i  %40s")
                     % i
                     % j
                     % pcps::format_number("%14.6f", _process_s_dot_s[i*_nsites+j] )
                  << std::endl;
        }
        outfile << std::endl;

        // write the S_i.S_j S_k.S_l averages
        outfile << boost::format("%6s  %6s  %6s  %6s  %40s") % "site i" % "site j" % "site k" % "site l" % "< ( S_i . S_j ) ( S_k . S_l ) >" << std::endl;
        for (int i = 0; i < _nsites; i++)
        for (int j = 0; j < _nsites; j++)
        for (int k = 0; k < _nsites; k++)
        for (int l = 0; l < _nsites; l++) {
          outfile << boost::format("%6i  %6i  %6i  %6i  %40s")
                     % i
                     % j
                     % k
                     % l
                     % pcps::format_number("%14.6f", _process_ssss[_nsites*(_nsites*(_nsites*i+j)+k)+l] )
                  << std::endl;
        }
        outfile << std::endl;

        // close the file
        outfile.close();

      }

  };

  //--------------------------------------------------------------------------------
  // pcps::PropertiesThreadAccum -- Accumulates correlation functions.
  //--------------------------------------------------------------------------------

  template <class S, class PART, class REF> class PropertiesThreadAccum : public pcps::ThreadAccum<S, PART, REF> {

    private:

      // pointer to the corresponding process accumulator
      pcps::PropertiesProcessAccum<S, PART, REF> * const _ppacc;

      // data members
      const int _nsites;            // number of sites in the lattice
      const int _configs_per_site;  // the number of configurations in a single site's Hilbert space
      S * _thread_sz;          // array to hold this thread's contribution to < Sz_i > for each site i
      S * _thread_sz_dot_sz;   // array to hold this thread's contribution to < Sz_i Sz_j > for each pair of sites i, j
      S * _thread_s_dot_s;     // array to hold this thread's contribution to < S_i.S_j > for each pair of sites i, j
      S * _thread_ssss;        // array to hold this thread's contribution to < S_i.S_j S_k.S_l>
      int * _workspace;             // array for use when computing local correlation contributions

    public:

      // constructor
      PropertiesThreadAccum(pcps::PropertiesProcessAccum<S, PART, REF> * const ppacc, const int tid, pcps::Arena & arena)
        : pcps::ThreadAccum<S, PART, REF>(ppacc, tid),
          _ppacc(ppacc),
          _nsites(ppacc->_nsites),
          _configs_per_site(ppacc->_configs_per_site)
      {

        // initialize arrays
        const S z = pcps::zero<S>();
        _thread_sz         = arena.allocate_array<  S>( size_t(_nsites)                         , z );
        _thread_sz_dot_sz  = arena.allocate_array<  S>( size_t(_nsites*_nsites)                 , z );
        _thread_s_dot_s    = arena.allocate_array<  S>( size_t(_nsites*_nsites)                 , z );
        _thread_ssss       = arena.allocate_array<  S>( size_t(_nsites*_nsites*_nsites*_nsites) , z );
        _workspace         = arena.allocate_array<     int>( size_t(_nsites)                         , 0 );

        // reset data
        this->reset();

      }

      // function to reset the accumulator
      void reset() {
        pcps::xscal(                        _nsites, pcps::zero<S>(), _thread_sz,        1);
        pcps::xscal(                _nsites*_nsites, pcps::zero<S>(), _thread_sz_dot_sz, 1);
        pcps::xscal(                _nsites*_nsites, pcps::zero<S>(), _thread_s_dot_s,   1);
        pcps::xscal(_nsites*_nsites*_nsites*_nsites, pcps::zero<S>(), _thread_ssss,      1);
      }

      // function to accumulate the correlation functions of the old configuration
      void accumulate(const long int iter, 
                      const long int count,
                      const bool accept,
                      const bool finished,
                      const S coeff_ratio,
                      const int * const new_config,
                      const int * const old_config,
                      pcps::CoeffValueFunc<S, PART, REF> & cvf,
                      pcps::ThreadAccum<S, PART, REF> ** accumulators,
                      const int n_accumulators) {

        if (accept || finished) {
          // compute correlation functions
          pcps::local_correlation_functions<S, PART, REF>(_nsites, _configs_per_site, count, old_config, _workspace, _ppacc->_avoRL, cvf,
                                                               _thread_sz, _thread_sz_dot_sz, _thread_s_dot_s, _thread_ssss);
        }

        if (finished) {

          // add this thread's contribution to the process's total
          pthread_mutex_lock(this->_pacc->mutex());
          const S total_samples = pcps::unity<S>() * this->_pacc->total_samples();
          pcps::xaxpy(                        _nsites, 1.0 / total_samples, _thread_sz,        1, &_ppacc->_process_sz.at(0),        1);
          pcps::xaxpy(                _nsites*_nsites, 1.0 / total_samples, _thread_sz_dot_sz, 1, &_ppacc->_process_sz_dot_sz.at(0), 1);
          pcps::xaxpy(                _nsites*_nsites, 1.0 / total_samples, _thread_s_dot_s,   1, &_ppacc->_process_s_dot_s.at(0),   1);
          pcps::xaxpy(_nsites*_nsites*_nsites*_nsites, 1.0 / total_samples, _thread_ssss,      1, &_ppacc->_process_ssss.at(0),      1);
          pthread_mutex_unlock(this->_pacc->mutex());

          // wait for other threads
          pthread_barrier_wait(this->_pacc->barrier());

          // combine results from all processors on root process
          if (this->_tid == 0) {

            const MPI::Comm & comm = MPI::COMM_WORLD;

            pcps::xcopy(                        _nsites, &_ppacc->_process_sz.at(0),        1, _thread_sz,        1);
            pcps::xcopy(                _nsites*_nsites, &_ppacc->_process_sz_dot_sz.at(0), 1, _thread_sz_dot_sz, 1);
            pcps::xcopy(                _nsites*_nsites, &_ppacc->_process_s_dot_s.at(0),   1, _thread_s_dot_s,   1);
            pcps::xcopy(_nsites*_nsites*_nsites*_nsites, &_ppacc->_process_ssss.at(0),      1, _thread_ssss,      1);

            pcps::xscal(                        _nsites, pcps::zero<S>(), &_ppacc->_process_sz.at(0),        1);
            pcps::xscal(                _nsites*_nsites, pcps::zero<S>(), &_ppacc->_process_sz_dot_sz.at(0), 1);
            pcps::xscal(                _nsites*_nsites, pcps::zero<S>(), &_ppacc->_process_s_dot_s.at(0),   1);
            pcps::xscal(_nsites*_nsites*_nsites*_nsites, pcps::zero<S>(), &_ppacc->_process_ssss.at(0),      1);

            pcps::reduce(comm, _thread_sz,        &_ppacc->_process_sz.at(0),                                _nsites, MPI::SUM, 0);
            pcps::reduce(comm, _thread_sz_dot_sz, &_ppacc->_process_sz_dot_sz.at(0),                 _nsites*_nsites, MPI::SUM, 0);
            pcps::reduce(comm, _thread_s_dot_s,   &_ppacc->_process_s_dot_s.at(0),                   _nsites*_nsites, MPI::SUM, 0);
            pcps::reduce(comm, _thread_ssss,      &_ppacc->_process_ssss.at(0),      _nsites*_nsites*_nsites*_nsites, MPI::SUM, 0);

            // compute diagonal terms of ssss correlations
            const double s = 0.5 * double(_configs_per_site-1);
            for (int i = 0; i < _nsites; i++)
            for (int k = 0; k < _nsites; k++)
            for (int l = 0; l < _nsites; l++) {
              _ppacc->_process_ssss[_nsites*(_nsites*(_nsites*i+i)+k)+l] = s * ( s + 1.0 ) * _ppacc->_process_s_dot_s[k*_nsites+l];
              _ppacc->_process_ssss[_nsites*(_nsites*(_nsites*k+l)+i)+i] = s * ( s + 1.0 ) * _ppacc->_process_s_dot_s[k*_nsites+l];
            }

          }

          // wait for thread 0
          pthread_barrier_wait(this->_pacc->barrier());

        }

      }

  };

  //------------------------------------------------------------------------------------------
  // pcps::PropertiesProcessAccum::create_thread_accum -- Function to create a thread accumulator
  //                                                      from a process accumulator.
  //------------------------------------------------------------------------------------------

  template <class S, class PART, class REF>
  pcps::ThreadAccum<S, PART, REF> * pcps::PropertiesProcessAccum<S, PART, REF>::create_thread_accum(const int tid, pcps::Arena & arena) {
    return new(arena) pcps::PropertiesThreadAccum<S, PART, REF>(this, tid, arena);
  }

} // end namespace pcps

#endif
